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1 Abstract 

We present a simple algorithm based on linear programming, that combines 
Kirchoff’s flux and potential laws and applies them to metabolic networks 
to predict thermodynamically feasible reaction fluxes. These law’s represent 
mass conservation and energy feasibility that are widely used in electrical 
circuit analysis. Formulating the Kirchoff’s potential law around a reaction 
loop in terms of the null space of the stiochiometric matrix leads to a simple 
representation of the law of entropy that can be readily incorporated into the 
traditional flux balance analysis without resorting to non-linear optimization. 
Our technique is new as it can easily check the fluxes got by applying flux 
balance analysis for thermodynamic feasibility and modify them if they are 
infeasible so that they satisfy the law of entropy. We illustrate our method by 
applying it to the network dealing with the central metabolism of Escherichia 
coli. Due to its simplicity this algorithm will be useful in studying large scale 
complex metabolic networks in the cell of different organisms. 
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Introduction 


Completely sequenced genomes open the door for wholistic understanding of 
metabolic networks. Study of metabolic networks is important in predicting 
the behaviour of the cell as a result of interactions between the different reac- 
tions occuring within the cell (Bailey, 1991; Palsson, 2000; Stephanopoulos 
et al. 1998). For example the flux balance analysis (FBA) has proved to 
be useful for studying the steady state metabolic flux inside the cell (Varma 
and Palsson, 1994). This approach is especially a valuable tool in the ab- 
sence of knowledge of detailed kinetic parameters of reactions inside the cell. 
Many authors have studied reaction networks from a very theoretical stand 
point (Oster et al. 1973; Peusner, 1986; Balabanian and Bickart, 1981), and 
also applied principles from non-equilibrium thermodynamics (Nicolis and 
Prigogine, 1977; Westerhoff and van Dam, 1987), but it has only been in 
the last few years, that these techniques have been applied to biological sys- 
tems (Bonarius et al., 1996, 1997; Pramanik and Keasling, 1997). Recently 
there have been attempts to further constrain the fluxes that are obtained by 
FBA, to satisfy the second law of thermodynamics by a method called energy 
balance analysis (EBA) (Beard et al. 2002). The EBA eliminates thermody- 
namically infeasible fluxes that are got by applying FBA alone (Price et al. 
2002 ). 

In this paper we incorporate the energy feasibilty constraint that is akin 
to Kirchoff’s potential law (Strang, 1986) into the linear FBA theory. This 
places additional constraints on the FBA solution space and eliminates ther- 
modynamically infeasible fluxes that do not satisfy the loop law. We illustrate 
our method by applying it to a part of the metabolic network of Escherichia 
coli (Delgado and Liao, 1997). This method will be useful in predicting the 
behaviour of large scale networks, especially in making predicitions of gene 
regulation and thermodynamic chemical potentials of the different chemical 
reactions that go on inside the cell. 

Previous attempts to formulate the EBA method resulted in a non-linear 
optimization problem (Beard et al. 2002), which can lead to errors if the 
method doesnot converge. The novelty of our approach is that we can solve 
it using linear programming, and that we can generate many solutions and 
test them for feasibilty. The approach taken in Beard et al. 2002 and Price et 
al. 2002 was to set fluxes to zero to satisfy the energy feasibilty constraints. 
This is a very restrictive approach as it misses out non-zero solutions. Our 
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method on the other hand first checks the fluxes without setting them to zero 
for feasibility and only sets fluxes to zero, if a sign transformation fails to 
make the flux distribution feasible. We also show mathematically that setting 
the flux of a reaction to zero, leads to no change in the chemical potential for 
that reaction. In Beard et al. 2002, nonzero change in the chemical potentials 
are computed for reactions with zero flux. They attribute reactions with zero 
fluxes to infinite resistance (or zero conductance). The objective of this paper 
is to present a simple technique that can be applied to large scale metabolic 
networks to predict the various internal fluxes and the corresponding change 
in the chemical potentials. To summarize we have developed a technique to 
test the thermodynamic feasibility of non-zero fluxes, and also to construc- 
tively generate feasible solutions. To accomplish this goal we develop some 
formal machinery using concepts from linear algebra. 


3 Flux Balance Analysis: Law of Flux Con- 
servation 

In a metabolic network inside the cell, several reactions catalyzed by different 
enzymes occur in concert. To study the rates of these reactions we resort to 
flux balance analysis. Applying the law of mass balance, the concentration 
of metabolites and the reaction rates (or fluxes) are related through the 
stoichiometry of the system. Using steady state approximation for the mass 
balance, the flux balance analysis has been formulated as a linear program 
by many authors (e.g. Varma and Palsson, 2000), in which there is a linear 
objective function, which could for example be to maximize growth, maximize 
ATP production, minimize glucose intake etc. This is written as a linear 
combination of the fluxes. 
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where, f E 71^ is the vector of n fluxes, S E 71^^^ is a stoichiometric matrix, 
m is the number of reactants in the network. All vectors by default will be 
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column vectors. Also, d, I and u are vectors £ 7^" of objective function 
coefficients, lower and upper bound constraints on the fluxes respectively, 
and 0 is a zero vector of size m. In equation (3) the inequalities are compo- 
nentwise for the vectors. The vector of objective function coefficients have 
to be determined experimentally, but in most cases has only one non-zero 
component corresponding to the flux one is trying to optimize. 

In the above formulation the number of fluxes n exceeds the number of 
metabolites m in the cell. So a convenient way to solve the above system of 
underdetermined equations is to resort to linear programming. Due, to the 
degeneracy in the problem, there are an infinite number of solutions possible 
that satisfy all of the constraints and optimize the chosen objective function. 


4 Energy Balance Analysis: Second law of 

Thermodynamics 

According to the second law, fluxes must flow from reactants of higher chem- 
ical potential to ones of lower chemical potential. This is the direction in 
which the entropy of the reaction increases (Qian et al. 2002). Since the 
FBA analysis gives rise to an infinite number of fluxes, many of these flux 
distributions violate the second law and hence are infeasible. From a network 
topology point of view, it is the presence of cycles in the flux direction, that 
violate the law of production of entropy. Applying Kirchoff’s loop law, one 
gets rid of these entropy violating cycles. 

From S we remove the columns corresponding to boundary fluxes and 
keep only the columns of non-redundant internal fluxes. The resulting matrix 
G G where, Uj is the number of internal fluxes in the network. Using 

the singular value decomposition (Strang, 1980) one can find the null space 
matrix N of G. The matrix N G consists of n/ orthonormal basis 

vectors of the null space of G. The dimension of the null space of G gives 
the number of independent and irreducible loops rii in the network (Strang, 
1986). By this we mean that a single basis loop cannot be decomposed 
into smaller loops. By taking linear combinations of these basis loops we 
can generate bigger and compund cycles. The singular value decomposition 
(SVD) is not unique for the G matrix, so the orthonormal basis for the null 
space of G is not unique. 
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Instead of obtaining the orthonormal basis for the null space of G from 
the SVD, one can obtain the rational basis for the null space from the reduced 
row echelon form of G. This basis is unique as the reduced echelon form of 
G is unique. 

Associated with each internal flux fi is a chemical potential difference 
A/Zi- These potential differences satisfy a law similar to the Kirchoff’s loop 
law in electrical circuits, namely (Beard et al. 2003), 


KAn = 0 (4) 

where, K = is a matrix whose rows are the basis vectors of the 

null space of G, and Ap G is a vector of chemical potential differences 
for the internal fluxes in the cell. 

The second law ensures that the entropy increases in each reaction i and 
hence the direction of flux fi is from metabolites of higher chemical potential 
to one of lower chemical potential and can be expressed as (Beard et al. 
2002). 


fAiJLi < 0 (5) 

Equations (4) and (5) are thermodynamic feasibility constraints that are 
applied in addition to the flux balance constraints. Equation (5) is a nonlinear 
constraint which when incorporated into the FBA makes the problem a non- 
linear programming problem. In this paper we propose a simpler algorithm 
to solve the problem as a linear programming problem. 

In addition to the above constraints one imposes upper and lower bound 
constraints on A^i. 


/3 < A^f < a (6) 

Where, /3 and a G represent the lower and upper bounds on the change 
in chemical potential A/^, and the inequality is componentwise. 

In the next section we describe the condition for checking the presence 
of cyclic fluxes, and we introduce some notation here. We will use upper- 
case indices to denote sets for example, let F be the set of all fluxes in the 
network, i? be the set of unrestricted fluxes, be the set of non-negative 
fluxes, and be the set of negative and positive fluxes respectively. 
Fluxes like /,• G F means that flux fi is a flux, is an unrestricted flux, ff^ is 
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a non-negative flux etc. The matrix N can be written in terms of its column 
vectors as N = where N*k is the A:th column vector 

of N. Also iV,fc = [nik,ri 2 k, -^riik, .■..,nmkf, where Uik is the (i,k) th entry 
of the matrix N. 


5 No cycle feasibility constraint 

In this section we introduce a simple test to detect the presence of loops 
in a metabolic network that violate the second law of thermodynamics. To 
do so we take advantage of the directionality of the flow of fluxes in the 
cycle. The number of rows rii of the K matrix gives the number of loops or 
cycles in the network (Strang, 1986). These loops are the basis cycles, as 
combining them gives bigger cycles. If in any row j of the K matrix Kj^ 
all the entries (there should be more than one non-zero entry) of the jth 
row are of the same sign, corresponding to the set of positive fluxes 
for the jth cycle, then the flux distribution is thermodynamically infeasible. 
If some unrestricted flux r, belonging to the cycle j is negative, it can be 
transformed to be positive by reversing the sign of the corresponding entries 
in the zth column of the G and K matrices, namely G*i and respectively. 
If after this transformation, any of the rows of the K matrix still have the 
same sign, then the flux distribution is thermodynamically infeasible, and 
we detect the presence of a cycle. By satisfying this condition we can get 
rid off the non-linear constraint in equation (5), and hence we transform the 
non-linear problem to a linear one. 

It should be noted that this no cycle feasibility condition is equivalent to 
solving the FBA problem with constraints (4) and (5). It can be seen that to 
satisfy equation (4) , for a single row of the K matrix corresponding to a single 
cycle, atleast one of the A/Zj should be of a different sign from the rest of 
the other components in the A/j, vector. This when combined with equation 
(5), prevents the formation of energy violating loops in the network, and 
hence satisfies the second law of thermodynamics. This no cycle feasibility 
condition is only applicable to basis loops, which can be easily observed in 
the echelon basis. 

Lemma 1: (sign transformation lemma) 

Transforming the unrestricted internal flux — fj to r* changes the sign of the 
ith column if*, of matrix K. 
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Proof: From the flux conservation equation Sf = 0, one can partition 
the flux vector / into internal flux vector x and boundary flux vector y. 
The columns of the stiochimetric matrix S can likewise be partitioned into 
columns G corresponding to the internal fluxes and columns H corresponding 
to the boundary fluxes. That is, 5 = [G' H]. The flux conservation equation 
can be rewritten as Gx = —Hy. 

Since N is the null space matrix of G, we have GN = 0, where 0 is a (m x nj) 
matrix of zeros. Consider the ith component of the unrestricted internal flux 
Tj, which is a component of x, that is negative. In the matrix vector product 
Gx, ri multiplies the ith column of matrix G. A negative r, value can 
be made positive by transferring the negative sign to all the elements in the 
column of G'*i. By this process the value of the matrix vector product Gx 
is unaffected. Hence G*i becomes — From the equation GN = 0 we 
have GN = where G*i, A'*! are the first 

columns of the G and N matrix respectively. This matrix-matrix product 
can be written compactly as where, there is implict sum- 

mation on the repeated index i. From this it is clear that transforming 
to — changes the sign of the ith row of matrix N, that is entries rin, ..., 
njn, change sign. Since by construction K = N^, the zth row of matrix N 
corresponds to the zth column K^i of matrix K, which then changes sign. 

By use of lemma 1 we transform all the unrestricted negative internal 
fluxes to reference positive fluxes, from which we can test the cycle condition 
for fluxes. 

Lemma 2; (zero transformation lemma) 

If the entries of the zth column of the matrix G, G*i, corresponding to the 
ith internal flux Xi, which only belongs to the ;th cycle are set equal to zero 
then the jth row of the K matrix Kj^ is zero everywhere except at the ith 
column, that is entry Kji is nonzero. Hence = 0. 

Proof: Since N is the null space matrix of G, we have GN = 0, where 0 is a 
(mxn;) matrix of zeros. From GN = 0, we have [G«inii....G*jny....G*inj„,] = 
0. Here there is summation on the repeated index i. The jth cycle corre- 
sponds to the jth row of K and hence the yth column of N matrix. Hence 
consider the equation corresponding to the jth column of N matrix 

G^iUij = 0 ( 7 ) 

where, there is summation over the repeated index i and 0 € is a column 
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vector of m zeros, and riij is (i, j)th element of the N matrix. 

G^lTlXj "H ... ~h G — G ^8} 

where on the left side of the above equation we exclude the ?th index as it is 
brought to the right. Without loss of generality let riij = 1, hence 

G '*!" ... “|~ G G (9) 
Since G*i is a zero vector, we have 

G*iTlij + ... + G^ni^nij = 0 ( 10 ) 

From the above equation (10) we see on the left hand side that, G,inij, ...., G*n; ’^no- 
are atmost TLi — 1 nonzero vectors, which were a part of the jth cycle along 
with GtiTiij (which corresponds to the ith internal flux), setting = 0 (is 
like setting Xi = 0) breaks the cycle, and the rest of the nonzero column 
vectors, G*i, ...., G*„. in the broken cycle are linearly independent by con- 
struction, since we only included non-redundant internal fluxes, and since 
atmost Tii internal fluxes are present in the jth cycle corresponding to the 
jth TOW of the if matrix Kj*, that are linearly dependent, breaking the cycle 
makes the rest of the fluxes in the jfth cycle linearly independent by construc- 
tion. The only way equation (9) holds is when nij = .... = n„.y = 0. Hence 
the jth colomn of the matrix N, is zero except for the riij element which 
is nonzero. Since K = the jth row of K, Kj^^ is zero except Kji = riij. 
Therefore, from equation (4) we see that Api = 0. 

From lemma 2 we see that when a particular internal flux is zero (we 
can zero out its respective column in the G matrix), then the corresponding 
change in chemical potential is zero, which is physically quite intuitive. This 
can be incorporated as an EBA constraint. Lemma 2 could be applied to 
overlapping cycles which share fluxes (two cycles are overlapping if they share 
atleast one flux, that is a column of the K matrix corresponding to the 
rational basis of the null space has two non-zero entries), in this case some 
cycles may be broken while some others may be created. It is therefore best 
to remove one cycle at a time by zeroing out fluxes that only belong to a 
single cycle, and are not shared among cycles. This is again got by observing 
the columns of the K matrix (corresponding to the rational basis of the null 
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space) that contain only a single nonzero entry corresponding to the flux that 
belongs to a single cycle. 

These lemmas provide an automated way of reducing the K matrix in 
large scale networks with many overlapping loops that share the same fluxes. 
The K matrix will also be called the feasibilty matrix. In the next section 
these concepts will be used in the EBA linear programming algorithm. 

6 EBA Algorithm 

1) Solve the FBA for the flux distribution. 

2) If all the fluxes are positive and the K matrix satisfies the no cycle feasibil- 
ity constraint, then the fluxes obtained from the FBA are thermodynamically 
feasible. We also find the corresponding A,u vector that satisfies constraints 
(4), (5) and (6). 

3) If some of the internal fluxes are negative then apply the sign transfor- 
mation lemma 1 to change the sign of the entries of the columns of the K 
matrix corresponding corresponding to these negative fluxes. Then if the K 
matrix satisfies the no cycle feasibility condition we get a feasible flux vector. 
The corresponding components of the A/u vector that satisfies constraints 
(4), (5) and (6) is computed. 

4) If after the above three steps, the K matrix is infeasible, we identify fluxes 
that are not shared by cycles. Corresponding to these fluxes, there is only 
one nonzero entry in their respective columns, which indicates the cycle in 
which the flux is present. For example we zero out any one flux Xi in some 
cycle j and apply the zero transformation lemma 2 to both the G and K 
matrix. 

In the K matrix then remove the jth row and ith column corresponding to 
the ;th cycle and ^th flux, and constrain A/ij = 0 and Xi = 0. Also remove 
the ith column from the G matrix. After this do steps (l)-(3) on the reduced 
system. If feasible compute A/j, vector that satisfies constraints (4), (5) and 
( 6 ). 

5) If not feasible step (4) can be repeated for each cycle one at a time and 
checking feasiblity until they are exchausted. 

In the above algorithm, computation of Ajj. is decoupled from the com- 
putaion of the flux distribution. In the MATLAB code however the two can 
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be combined into an equivalent linear programming formulation, the details 
are not discussed in this paper. Also, when trying to identify non-shared 
fluxes, it is best to use the rational basis of the null space for the K matrix. 

7 Application of EBA Algorithm 

In this section we consider the example discussed in Beard et al (2002). The 
set of reactions are shown below. 

rxn 1; A + 2B < > C 

rxn 2: C + D < — > 2A + 2B 
rxn 3: A + B < — > 2D 

rxn 4: A + C < > B 4- 3D 

rxn 5: B < — > D 

This reaction network has 5 internal fluxes corresponding to the 5 reac- 
tions, and associated with each internal flux is a change in chemical potential. 
In addition to these there are 4 boundary or exchange fluxes corresponding 
to the 4 metabolites A, B, C and D. 

For this reaction network we wish to determine the maximum steady- 
state production of reactant D, for a given maximal input flux of reactant C. 
This problem assumes that reactant C is the only available input substrate, 
and that its value is set equal to 1, and that for A and B their values are set 
equal to 0. 

The stoichiometric matrix G for the four metabolites A, B, C and D cor- 
responding to the four rows respectively and the flve reactions corresponding 

0 ■ 

-1 
0 

1 . 

four metabolites as rows and 
four metabolites respectively 

as the four columns is 


to the flve columns of G respectively, is 


G = 


-1 2 -1 

-2 2 -1 

1 -1 0 

0-1 2 


-1 

1 

-1 

3 


Also, the stoichiometric matrix H for the 
the four exchange fluxes corresponding to the 
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10 0 O' 

0 10 0 

0 0 1 0 

0 0 0 -1 . 

The full stoichiometric matrix is 5 = [(? H] and the flux distribution 
/ 6 7^® that maximizes the production of D, and that satisfies the flux 
conservation constraint Sf = 0. The five internal fluxes corresponding to 
the five internal reactions and the boundary flux corresponding to D are 
unrestricted, while the other three boundary fluxes corresponding to A, B 
are less than or equal to 0, and that for C it is less than or equal to 1. 

The orthonormal null space vectors corresponding to matrix G are given 
in Beard et al. (2002) and can be computed easily using MATLAB, and are 

kx^ = [-0.7163, -0.3205, 0.4710, -0.3958, -0.0752] (11) 

= [-0.3345, -0.4347, -0.6349, 0.1001, 0.5348] (12) 

The dimension of the null space of matrix is 2. The K matrix consists 

nr nr 

of kx and ^2 as its rows. 

Alternatively, one can use the rational null space vectors corresponding 
to matrix G, that are 

fcj = [2 ,l,-l,l, 0] (13) 

fc2j = [-1,-1, -1,0,1] (14) 

The K matrix then consists of and k 2 ^ as its rows. This representation 
tends to give a clear picture about shared and non-shared fluxes among cycles. 
Here fluxes Xi , X 2 and 0:3 are shared among the two cycles, whereas X 4 belongs 
to cycle 1 and X5 belongs to cycle 2. 

Optimal flux vector f = [x^, 6 and change in chemical potential 

vector Ay, e are obtained using the standard linear programming routine 
in MATLAB. They are 


H = 


y = [0,0,l,3f ^5) 

X = [-0.154, 0.445, 0.643, 0.401, 0.956]^ (16) 

Ay = [0.6178, -0.8158, -0.5801, -1.000, -0.7780f (17) 
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Where x € 7^^ is a flux vector of the 5 internal fluxes and y € 72.“^ is a 
flux vector of the 4 boundary fluxes. From equation (15) we observed that 
3 units of D are produced per unit of C consumed. All the fluxes have been 
normalized with respect to one unit of C. 

The optimal flux distribution obtained using our algorithm that satisfies 
the thermodynamic constraint is the same as the one obtained by Beard et 
al (2002). But according to Beard et al. it is infeasible thermodynamically. 
But according to our sign test we see that it is feasible and we also provide a 
corresponding vector for the x internal flux distribution, that satisfies all 
the feasibility constraints given in equations (4) and (5). Since in this exam- 
ple has unrestricted components, in constraint (6) the lower and upper 
bounds are -oo and -l-oo respectively. The formulation given in Beard et al. 
might have failed to converge to a A^ vector and hence the flux distribution 
was reported infeasible. Our criterian for the EBA test is quite robust and we 
can immediately tell if the flux distribution is thermodynamically infeasible, 
without carrying out any nonlinear optimization. 


8 Analysis of F. Coli Central Metabolism 

We use the stoichiometric matrix S of the model E. coli system from Table 1 
(Delgado and Liao (1997)) for our FBA/EBA analysis. The reaction network 
contains 19 metabolites linked by 23 reactions (Figure 1 Delgado and Liao 
(1997). Out of these 23 fluxes there are 3 external or boundary fluxes and 
the rest 20 are internal fluxes. The network considered takes glucose as input 
and produces acetate and carbon dioxide. The energy and the metabolites 
involved in this process are used for the synthesis of proteins, DNA, RNA etc. 
We applied our algorithm to maximize the production of biomass flux, which 
is a linear combination of the different fluxes with experimentally determined 
stoichiometric coefiicients, and are given in Neidhardt et al. (1990). These 
coefficients are for the conversion of key metabolites to biomass of E. coli 
propotoplasm. In the FBA optmization the internal fluxes are unrestricted, 
and only satisfy the flux balance constraint. The C02 and acetate fluxes 
come from the literature (Kleman and Strohl (1994); Shiloach et al. (1996)). 
Since only the relative rates matter, the glucose flux is set to 1, and all other 
fluxes are normalized with respect to it. 

The G matrix is formed by considering the columns of the following in- 


12 



ternal fluxes from Table 1 of Delgado and Liao (1997); 

^ ~ \Jpgii Jpepi Jpyki Jpdhi Jacet Jzi Jicti •Al) Jl2 ) Jppcj Jl4i Jl5j Jl6t Jtkti 
Jtai, Jresp, Jatp, Jbiomass, Jgiyox], and the H matrix is formed from the columns 

of the external fluxes j/ = [Jjiucj^cozt^ace]- 

The null space of the G matrix is of dimension 1, hence the K matrix consists 
of one row. This also indicates the presence of a single loop in the network. 

K = [0, 0, 0, 0.24, 0.24, 0, 0, -0.24, -0.25, -0.24, -0.25, 0, 0, 0, 0, 0, 

— 0.24, —0.73, 0, 0.24]. We can also use the rational null space vector k\r of 

G to construct the matrix K as before 

fc J = [0, 0, 0, 1, 1, 0, 0, -1,-1, -1, -1, 0, 0, 0, 0, 0, -1, -3, 0, 1]. 

From the above equation we see the following 9 fluxes which form a cycle: 
\Jpyki Jpdhj Jicti J\2i Jppci Jrespi Jatpt Jgiyox\- These Correspond to the non- 
zero entries in the row of the K matrix and also in the kx^ vector. 

The optimal flux vector / = [x^,y^Y ^ change in chemical 

potential vector satisfying the flux balance and thermodynamic 

constraints computed by our algorithm are: 

2/ = [l,2.2,0.3f 

X = [0.87, 0.85, 1.58, 1.92, 2.71, 0.27, 0.56, -1.03, -1.11, -1.11, -1.38, 

0.12, 0.09, 0.03, 0.03, 0.03, 1.80, 2.95, 0.0001, 1.59f 

Aa/-[- 8.25, 8.25, 8.25, 6.27,-6.27,-0.25,-8.25,10.95,10.95,10.95, 
10.95, -8.25, -8.25, -8.25, -8.25, -8.25, -10.74, -17.30, -8.25, -6.27]^ 
Where x 6 is a flux vector of the 20 fluxes including the 19 internal 
fluxes and the biomass flux, and y € is a. flux vector of the 3 boundary 
fluxes. In this example we have rounded the numerical values to two decimal 
places. The optimized biomass flux is Jbiomass = 7-27 x 10“® « 0.0001 per 
unit of glucose consumed. 

To see if the x vector of internal fluxes is thermodynamically feasible, 
we apply the sign transformation lemma to the negative flux components 
xs,Xg,XiQ and Xn, and see that the K matrix satisfies the feasibility criteria 
discussed in the previous sections. 


9 Conclusions 

In this paper we give a simple linear programming algorithm for the flux and 
energy balance analysis. Instead of applying a nonlinear thermodynamic fea- 
sibility constraint, as has been done previously, it uses the sign of the null 
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space to decide if the flux distribution of the metabolic network computed 
by flux balance anaysis satisfies the second law of thermodynamics. This 
technique is different from the previous approaches, as it is constructive, and 
can generate several feasible solutions for the metabolic network. A self con- 
sistent connection between setting a reaction flux to zero with respect to a 
corresponding change in the chemical potential for that reaction is developed. 
This can be used as a energy feasibility constraint for the reaction. We ap- 
plied the method to a part of the metabolic network of E. coli and computed 
the fluxes and change in chemical potentials for the internal reactions. It 
should however be noted that FBA together with EBA are still not able to 
constrain the metabolic network completely and this leads to an infinity of 
flux and change in chemical potential distributions. More realistic bounds 
on the values of fluxes and change in chemical potentials are required to fur- 
ther constrain the system. This could be got from studying the biochemistry 
of several pathways. The connection between linear algebra and thermody- 
namics is exploited to formulate the law of entropy in metabolic networks. 
Furthermore, as the analysis and technique presented here are simple and 
straightforward, they can easily be applied to study large scale metabolic 
networks. 
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